{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Vortex Sheet Roll-up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "using PotentialFlow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "using Plots\n",
    "clibrary(:colorbrewer)\n",
    "default(legend = :none, ratio = 1, linewidth = 2,\n",
    "        grid = false, colorbar_title = \"\\$\\\\gamma\\$\", cbar = :right,\n",
    "        clim = (-0.5, 0.5), xlim = (-1.2, 1.2));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Initializing the sheet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We want to recreate the simulation of a vortex sheet roll-up as presented in [Krasny 1987: Computation of vortex sheet roll-up in the Trefftz plane](https://doi.org/10.1017/S0022112087002830).\n",
    "The discretization parameters are listed below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "N = 200\n",
    "α = linspace(0, π, 2N+1)\n",
    "zs = complex.(cos.(α))\n",
    "Γs = sin.(α)\n",
    "δ = 0.05;\n",
    "Δt = 0.01;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We can initialize the vortex sheet by providing it with a list of control points `zs`, the circulation from the tip of the sheet `zs[1]` to any subsequent control point, and a blob radius.\n",
    "\n",
    "In order to use the time marching schemes defined in `src/TimeMarching.jl`, we define `compute_ẋ!` that computes the velocities of the sheet control points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "compute_ẋ! (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function compute_ẋ!(ẋ, x, t)\n",
    "    reset_velocity!(ẋ, x)\n",
    "    self_induce_velocity!(ẋ, x, t)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Time Marching"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We start by using forward-Euler:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "sheet = Vortex.Sheet(zs, Γs, δ)\n",
    "w = allocate_velocity(sheet);\n",
    "\n",
    "for t in 0.0:Δt:1.0\n",
    "    forward_euler!(sheet, sheet, t, Δt, compute_ẋ!, advect!, w)\n",
    "    # which does the same thing as:\n",
    "    #     reset_velocity!(w, sheet)\n",
    "    #     self_induce_velocity!(w, sheet)\n",
    "    #     advect!(sheet, sheet, w, Δt)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(sheet, color = :RdBu, ylim = (-0.5, 0.2), size = (600, 200))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "which seems to work just fine.  But over a longer time interval, instabilities in the sheet start to show, resulting in sheet segments that overlap with each other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "for t in 1.01:Δt:4.0\n",
    "    forward_euler!(sheet, sheet, t, Δt, compute_ẋ!, advect!, w)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(sheet, color = :RdBu, ylim = (-1.4, 0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "If we use 4th-order Runge-Kutta, as was done in [Krasny 1987](https://doi.org/10.1017/S0022112087002830), we obtain:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "sheet  = Vortex.Sheet(zs, Γs, δ)\n",
    "sheet₊ = deepcopy(sheet)\n",
    "\n",
    "ẋs = [allocate_velocity(sheet) for k in 1:4]\n",
    "for t in 0:Δt:4.0\n",
    "    rk4!(sheet₊, sheet, t, Δt, compute_ẋ!, advect!, ẋs)\n",
    "    sheet₊, sheet = sheet, sheet₊\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(sheet, color = :RdBu, ylim = (-1.4, 0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "The higher order time-marching scheme gives us a more stable vortex core, but the large strain-rates start to reveal the discrete nature of our vortex sheet.\n",
    "\n",
    "Another way we can supress the instabilities is to use discrete Fourier transforms to directly filter out any length scales in the sheet below a specified threshold.  As an added benefit, the filtering process automatically inserts and redistributes control points to maintain a uniform particle density along the sheet.\n",
    "\n",
    "Here, we will filter out any distrubances with length scales greater than 0.4 times the blob radius and redistribute the control points so that there are 5 control points within a blob radius.  If we use forward-Euler with filtering, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "sheet₊ = Vortex.Sheet(zs, Γs, δ)\n",
    "ẋ = allocate_velocity(sheet)\n",
    "for t in 0:Δt:4.0\n",
    "    forward_euler!(sheet₊, sheet₊, t, Δt, compute_ẋ!, advect!, ẋ)\n",
    "    z₌, Γ₌, l = Vortex.Sheets.remesh(sheet₊, 0.2δ)\n",
    "    Vortex.Sheets.filter_position!(z₌, 0.4δ, l)\n",
    "    Vortex.Sheets.redistribute_points!(sheet₊, z₌, Γ₌)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(sheet₊, color = :RdBu, ylim = (-1.4, 0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "and if we use 4th-order Runge-Kutta on top of the filtering, we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "sheet₋ = Vortex.Sheet(zs, Γs, δ)\n",
    "sheet₊ = deepcopy(sheet₋)\n",
    "\n",
    "ẋs = [allocate_velocity(sheet₋) for k in 1:4]\n",
    "for t in 0:Δt:4.0\n",
    "    rk4!(sheet₊, sheet₋, t, Δt, compute_ẋ!, advect!, ẋs)\n",
    "    Vortex.Sheets.filter!(sheet₊, 0.2δ, 0.4δ)\n",
    "    sheet₊, sheet₋ = sheet₋, sheet₊\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(sheet₊, color = :RdBu, ylim = (-1.4, 0.2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "The continuity of the sheet comes at the cost of increased particle count.  Before filtering, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "401"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "length(sheet.blobs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "with filtering, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1661"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "length(sheet₊.blobs)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1-pre",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
